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Abstract 

In recent computer simulations of a simple monatomic system interacting 
via the Dzugutov pair potential, freezing of the fluid into an equilibrium 
dodecagonal quasicrystal has been reported [M. Dzugutov, Phys. Rev. Lett. 
70, 2924 (1993)]. Here, using a combination of molecular dynamics simulation 
and thermodynamic perturbation theory, we conduct a detailed analysis of the 
relative stabilities of solid-phase structures of the Dzugutov-potential system. 
At low pressures, the most stable structure is found to be a bcc crystal, which 
gives way at higher pressures to an fee crystal. Although a dodecagonal 
quasicrystal and a cr-phase crystal compete with the bcc crystal for stability, 
they remain always metastable. 

PACS numbers: 61.43.Bn, 61.44.Br, 61.50.Ah, 64.70.Dv, 64.70.Pf 
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I. INTRODUCTION 



Several years ago [|TJ, Dzugutov introduced a new model pair potential for the purpose 
of studying glass transitions by means of molecular dynamics (MD) simulation. Most sim- 
ulations of glass transitions to date have been performed using binary mixtures, since one- 
component simple liquids, when supercooled, readily nucleate crystallites. The Dzugutov 
potential, which features a pronounced maximum at a range typical of next-nearest-neighbor 
coordination distances in close-packed crystals, suppresses crystallization by construction, 
and thus facilitates glass formation. Following its initial successful use in studies of su- 
percooled liquids |],|2|], the Dzugutov potential was subsequently adopted in simulations of 
freezing HQ]. Contrary to expectations, however, the observed solid structure was deter- 
mined to be, remarkably, a monatomic dodecagonal quasicrystal. The structure, also known 
as tetrahedrally or topologically close packed (tcp) ||, is of the Frank-Kasper type and is 
composed of layers of square-triangle tilings. In previous work, the Dzugutov potential and 
the associated dodecagonal structure have been used to study self-diffusion in quasicrys- 
tals §. 

The motivation for the present work stems primarily from our interest in the nucleation 
and stability of quasicrystal phases. Further motivation comes from the realm of colloid 
physics, where, given the extreme tunability of colloidal interparticle interactions, it is con- 
ceivable that a Dzugutov-like pair potential might be engineered to produce bulk samples of 
one-component quasicrystals fTj. The main purpose of the study reported here is to chart, 
by means of both MD simulation and thermodynamic perturbation theory, the fluid-solid 
phase diagram of the Dzugutov-potential system. 

The paper is organized as follows. In Section |J, after specifying the pair potential, 
we give details of the simulation methods, and describe our analyses of the resulting solid 



structures. Section |T| outlines the theoretical methods used, while Section [TV] characterizes 
the quasicrystal and other structures of interest. In Section |V| we present results - from both 
simulation and theory - for relative stabilities of competing solid structures. Among our 



main results, we find that (1) at low temperatures and pressures, the stable solid structure 
is the bcc crystal; (2) at high pressures, the stable solid is the fee crystal; and (3) the stable 
ground state (T = 0) is never the dodecagonal quasicrystal, but rather a periodic crystal 
whose structure depends on the pressure. Finally, in Section [VI] we summarize and discuss 
implications of the results for future work. 



II. MOLECULAR DYNAMICS SIMULATIONS 



A. The Interaction 



The Dzugutov pair potential 0], plotted in Fig. 1, is defined by 
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The potential is characterized by a minimum at r = 1.13 a of depth —0.581 e, having 
the same form as that of the Lennard- Jones potential, followed by a maximum at r = 1.63 
a of height 0.460 e. The maximum is designed to prevent the system from crystallizing into 
simple crystal structures. Beyond the maximum the potential tends to zero continuously and 
is cut off at a range of r c = 1.94 a, which ensures that CPU times remain within reasonable 
limits. 



B. Simulation Method 



Classical isothermal (NVT) and isothermal-isobaric (NPT) MD simulations were per- 
formed using the constraint method || . An extension of this method allows us to introduce 
constant temperature or pressure gradients ||. Newton's equations of motion were inte- 
grated using a fourth-order Gear predictor-corrector algorithm (see, e.g., ||) with a time 
increment of St = 0.0005 a^Jm/e for all simulations. Periodic boundary conditions were 
applied to an orthorhombic simulation cell, whose volume (in the NPT simulations) was 
permitted to change isotropically. 

The sample sizes range from 54 to 1024 atoms, with the lengths of the cell along the 
three orthogonal coordinate axes being chosen to make the sample shapes as close to cubic 
as possible. Simulations with the stable phases were carried out with samples containing 54, 
250, and 1024 atoms for bcc; 60 and 480 atoms for the cx-phase; and 108 and 500 atoms for 
fee. Most results are reported for 250, 500, and 1024 atoms, though if not explicitly stated, 
the sample size is irrelevant. The potential energy and enthalpy per atom were found to be 
quite insensitive to sample size. 

C. Structural Analysis 

To analyze the structures that arise in a simulation at finite temperature, we quenched 
the system by setting the temperature in our NPT-MD program to zero, thereby using the 
program as a steepest descent algorithm. Quenching forces the system to seek out its local 
energy minimum, which facilitates structure identification. 

Dzugutov || has pointed out that vacancies may play a major role in stabilizing qua- 
sicrystalline structures. To determine the number of vacancies in a sample, we begin by 
constructing the Voronoi cells and their dual, the Delaunay cells, and determining from the 
latter the distribution of free volumes in the structure. Although the free volume distribu- 
tion gives a reasonable representation of the interstitial sites, it greatly overestimates the 



vacancies (by a factor of about ten). This is because the Delaunay cells are face-to-face 
packed tetrahedra, whereas the vacancies should be represented by spheres that could easily 
cover several tetrahedra. To determine the correct number of vacancies, we first select the 
Delaunay cells that are large enough to accommodate an atom and do not overlap with an 
already existing atom. With these cells we create trees of mutually overlapping cells, the 
nodes of the tree being the centers of the Delaunay cells and the edges the distance vectors 
between cells too close to be filled simultaneously. We then fill the tree with spheres, starting 
at the outermost ends. After adding a sphere, all the Delaunay cells connected to its node 
are discarded. The next sphere is then added onto the next outermost node remaining. The 
procedure is repeated until the whole tree is filled, after which the algorithm repeats with 
another Delauney cell not belonging to the current tree. This method allows us to fill the 
sample as densely as possible with vacancies. 

In order to help characterize and distinguish the solid structures observed in the simula- 
tions, we have computed from the atomic coordinates both radial and angular distribution 
functions. However, such averaged functions often do not allow unique identification of a 
structure, which requires as well the spatial distribution of bonds. Thus, we have also gen- 
erated bond order diagrams. These are stereographic projections of the nearest-neighbor 
bonds constructed as follows. First, all neighbor vectors are determined and normalized to 
unit length. Next, the vectors are placed at a common origin so that their endpoints lie 
on the unit sphere. Finally, the distribution of the points on the sphere is represented by 
stereographic projections along the three coordinate axes. The pictures thus obtained reveal 
the global symmetry of the sample. 

III. THEORY 

For comparison with the simulation data, we have independently calculated the phase be- 
havior of the system by means of thermodynamic perturbation theory. Taking the Dzugutov 
pair potential as input, we apply the approximate theory of Weeks, Chandler, and Andersen 



(WCA) [|TIJ to a classical system of iV pair-wise-interacting particles in a volume V. The 
WCA approach is especially well suited to pair potentials that contain a steeply repulsive 
core and has been successfully applied to the Lennard- Jones potential |TT|-p^], which has a 
repulsive core similar to that of the Dzugutov potential. 

The WCA approximation splits the pair potential 0(r) at its first minimum into a short- 
range repulsive reference potential (po(r) and a perturbation potential 4> p {r) and prescribes a 
mapping of the reference system onto an effective hard-sphere (HS) system. The Helmholtz 
free energy F of the system separates correspondingly into reference and perturbation parts, 
according to 

F = F + f 1 dA($ p )A, (4) 
J o 

where Fq is the free energy of the reference system, 

$ p = ^0 p (| ri -r,|) (5) 

i<j 

is the total perturbation energy, and (• • -)\ denotes averaging with respect to the probability 
distribution of a system with pair potential <f>\(r) = <po{r) + \<fi p (r). Expansion of ($ p )a 
in powers of A about the reference system (A = 0) generates an exact perturbation series. 
Mapping the reference system onto an effective HS system, the free energy may be expressed, 
to first-order in the perturbation potential, as 

F[p(r)] = F HS [p(r)] + — — / dr' r' 2 g m (r'; [p(r)]) </> p (r'), (6) 

V Jo 

where Fns[p(r)] and gHs( r ; [p( r )D are the f ree energy and radial distribution function (rdf), 
respectively, of the HS reference system, both functionals of the equilibrium one-particle 
number density p(v). The rdf is defined, in turn, according to 

g us (r; [p(r)]) = J dQ f dr'p (2 )(r', r' + r), (7) 

as an orientational and translational average of the two-particle density p*- 2 ^(r, r'). The 
second- and higher-order terms are proportional to successively higher powers of inverse 



temperature 1/T, the coefficients being related to mean fluctuations of $ p [Q. Accuracy of 
the first-order approximation [Eq. @] is thus assured as long as fluctuations in $ p remain 
sufficiently small relative to the thermal energy k&T. 

The free energy of the fluid phase is calculated via the uniform limit (p(r) — > p) of 
Eq. (El), using the essentially exact Carnahan-Starling and Verlet-Weis forms jnj for the 



HS free energy per particle, fns(p), and rdf, gns( r ), respectively. For the solid phase, the 
HS free energy functional is approximated by means of classical density-functional (DF) 
theory jn|. The DF approach is based on the existence of a functional jF[p(r)] of the 
density p(r) that satisfies a variational principle, according to which .F[p(r)] is minimized - 
for given average density and external potential - by the equilibrium density, its minimum 
value equaling the Helmholtz free energy F. In the absence of an external potential, JF[p(r)] 
may be decomposed into an exactly known ideal-gas contribution 

^id[p(r)] = k B T J drp(r) [ln(p(r)A 3 ) - l] , (8) 

which is the free energy in the absence of interactions (A being the thermal de Broglie wave- 
length) and an excess contribution jF ex [p(r)], depending entirely upon internal interactions. 

Here we approximate the excess free energy of the HS solid by the modified weighted- 
density approximation (MWDA) H16J17 1, which gives a reasonable description of the HS 



system. The MWDA maps the excess free energy per particle of the solid onto that of a 
corresponding uniform fluid of effective density, according to 

L^MWDA W)] = /Hg(/3); (g) 

where the effective (or weighted) density 

P = JjJ dr J dr'p(r)p(r')w(\r-r'\;p) (10) 

is a self-consistently determined weighted average of p(r). The weight function w(r) is 
specified by normalization and by the requirement that ^F^ WDA [pij)] generate the exact 
two-particle (Ornstein-Zernike) direct correlation function c(r) in the uniform limit. This 



leads to an analytic relation [16| between w(r) and the fluid functions fns and c(r), computed 
here using the solution of the Percus-Yevick (PY) integral equation for hard spheres 



Practical calculation of Fns[p(r)] an< l 9ns( r ', [p( r )D requires specifying the solid density, 
i.e., the coordinates of the lattice sites (equilibrium particle positions) and the shape of 
the density distribution about these sites. Here we consider fee, hep, bec, and a-phase 



crystals and the dodecagonal quasicrystal structures described below in Sec. |IV|- The density 
distribution is modelled by the Gaussian ansatz. This places at each site R a normalized 
isotropic Gaussian, such that 

p( r ) = (^) 3/2 E ex p(-«l r -RI 2 )> (ii) 

R, 

the single parameter a determining the width of the distribution. The Gaussian ansatz 
has been shown by simulation [[18] to reasonably describe the density distribution of close- 
packed crystals. For nonoverlapping neighboring Gaussians - consistently the case here - 
the ideal-gas free energy per particle [Eq. (H)] is very accurately approximated by 

^id = ^ B Tm(aA 2 )-| (12) 

to within an irrelevant constant. The HS free energy is obtained, for a given solid structure 
and average density, by minimizing the approximate functional J-Hs[p( r )] = -7"id[p(r)] + 
F™ WDA W)\ [ fr o m E qs- ©, (23), and (|TJ)] with respect to a. Predictions of the MWDA 
for free energies and pressures of HS solids are in good agreement with simulation data for 
both fee [16] and bec [[HJ crystals. Although the theory underpredicts, by roughly 20 %, 
the Lindemann ratios (ratio of root-mean-square particle displacement to nearest-neighbor 
distance) at melting for both crystal symmetries, it is only the free energies that determine 
thermodynamic phase behavior. Note that in simulations of the HS bec crystal a constraint 
of single-cell occupancy is usually imposed in order to stabilize the crystal against shear. 

The perturbation free energy in Eq. @ requires knowledge of the hard-sphere rdf, which 
may be expressed, in general, as a sum over coordination shells: 

oo 

Wr;[p(r)]) = £s W (r;[p(r)]). (13) 

i=l 



The functional g^(r; [p(r)]) are obtained using the approach of Rascon et al. [I3|, which 



approximates the second and higher coordination shells in mean-field fashion and corrects 
the first coordination shell for nearest-neighbor correlations. Thus, ignoring correlations for 
i > 2, and substituting p^(r, r') = p(r)p(r'), together with Eq. (|TT[) , into Eq. (^), yields 

2»(r; [p]) = ^(|^) 1/2 ^" ex P[-«^ - RifIA * > 2, (14) 

where rii is the coordination number and Ri the lattice vector magnitude of the ith shell. 
The first peak is parametrized by 

m/ r / m\ Aexp[-ai(r - ri) 2 /2] 
#«(r;[p(r)]) = ^ £ >-LA, r > d, (15) 

where d is the effective HS diameter and where the parameters A, a\, and r\ are determined 
by imposing three sum rules, namely the virial equation (relating the contact value to 
the bulk pressure P), normalization to the nearest-neighbor coordination number rii, and 
approximation of the first moment by its mean-field value. Together then, the HS pressure 
P = p 2 dfus/dp and the value of a that minimizes Fns[p(r)] determine <7hs(?"; [p( r )D an d 
so the perturbation free energy for a given solid structure. The approximation expressed 
by Eqs. (0)-(EO!) ^ s m excellent agreement with simulation data for the HS fee crystal, and 
has been successfully applied, in a perturbation theory, to Lennard- Jones and square-well 
solids The approximation also has been tested against, and found to closely match, 

Monte Carlo simulation data for ^Hs( r ; [p( r )D of a HS bec crystal (at density pa 3 = 1.1), 
subjected to a single-cell occupancy constraint to suppress shear instability JZ(J. Further 
simulations will be required, however, to test the approximation for the HS bec crystal 
at lower densities, where next-nearest-neighbor correlations and anisotropies in the density 
distribution may not be negligible. 

It remains still to specify the effective HS diameter d. According to the WCA prescrip- 
tion, d is the root of the nonlinear equation 

J dr WHS (r;[p(r)];d)Ae(r) = 0, (16) 



where yns{r] [p(r)];d) = exp[0ns( r ; ^)/^B^]fi'Hs( r ; [p( r )]id) is the HS cavity function and 

Ae(r) = exp[-Mr)/k B T] - exp[-0 HS (r; d)/k B T] (17) 

is a function that is nonzero only over a narrow range £d (£ <C 1) around r = d. This 
choice ensures that the free energy of the reference system differs from that of the effective 
HS system only by terms of 0(£ 4 ) and higher. In practice, lacking knowledge of the cavity 
function of the HS solid for r < d, we expand the quantity r 2 y H s( r i [p( r )];^) i n a Taylor 
series about r = d and retain the first three terms. 

The theory set out above provides a reasonable approximation for the Helmholtz free 
energy of the system at temperatures of order e/k B and higher. From the free energy, any 
bulk thermodynamic property then may be calculated. Of particular relevance to phase 
behavior are the pressure and chemical potential. In Sec. [VD| , we present our theoretical 
predictions for the phase diagram of the Dzugutov-potential system. 

IV. SOLID STRUCTURES 

A. Dodecagonal Quasicrystal 

The structural model of the dodecagonal quasicrystal that we have investigated is a 
layered system that is periodic in one direction, but quasiperiodic and twelve-fold symmetric 
in the perpendicular plane |2l[| . It is of Frank-Kasper type, i.e., it is mostly tetrahedrally 
close-packed, and can be described as a periodic ABAB stacking of a primary dodecagonal 
layer A and two secondary hexagonal layers, B and B, which are rotated by 30° with respect 
to each other to obtain dodecagonal symmetry. The atoms in layer A form the vertices of 
a simple tiling made of squares, triangles, 30° rhombi and two kinds of hexagons. The 
threefold symmetric hexagon is known as the "shield". These tiles, together with their 
decorations, are shown in Fig. |2|. A sample of a square-triangle configuration is displayed in 
Fig. |3|. The dodecagonal quasicrystal structure also can be regarded as the decoration of a 
simple dodecagonal tiling p2|f23"| . 



The stability of the monatomic Frank-Kasper-type decoration of the square-triangle- 
rhombi-shield tiling with the potential of Eq. ([!]) was reported by Dzugutov ||. Upon 
cooling below the glass transition temperature, a glass forms, which transforms, after a very 
long annealing time, into a dodecagonal quasicrystal. The underlying tiling structure is 
mainly a decorated square-triangle tiling with a few rhombi and shields. 

The aperiodicity of quasicrystals forbids periodic boundary conditions in the simulation. 
Taking a finite patch with open boundary conditions also should be avoided, as the surface 
energy would affect structural stability. A solution is to use periodic approximants, which are 
finite, orthorhombic cells whose boundaries match on opposite sides. In this way, periodic 
boundary conditions may be applied, as is done throughout this paper. 

B. Square-Triangle Crystals 

In addition to the quasicrystalline tilings, it is also possible to generate crystalline phases 
with squares and triangles decorated in the same fashion as for the quasicrystals. If only 
squares are used, the A15- or /5-W-phase (also known as cP8 Cr 3 Si) is obtained, whereas a 
pure triangle tiling results in the Z structure (hp7 Zr 3 Al 4 ). If both squares and triangles are 
permitted, the cx-phase or /3-U (also known as tp30 Cr 46 Fe54) and the H-phase are obtained. 
The two phases differ in the arrangement of the tiles. The vertex configurations of the 
crystalline phases are shown in Fig. |j. They are denoted by a, z, h, a according to the phase 
in which they appear. 

The unit cell of the a-phase can be subdivided into two regular triangles and two squares. 
The atoms at the vertices of these tiles are 14-fold coordinated, while the atom in the center of 
the triangles is 15-fold coordinated. The remaining atoms on the edges and in the interior of 
the square are 12-fold coordinated icosahedra. Since all coordination shells have a triangular 
surface, all atoms are tetrahedrally close packed. Pure square-triangle structures and tilings 
with additional shields are very stable Rhombi, however, are unstable and transform 
into the other tiles. 



V. RESULTS 



A. Analysis and Comparison of Solid Structures 

In this section we discuss the geometric properties of the amorphous structure, the bcc 
and fee crystals, the nucleated tcp-phase, and the cx-phase. We have used three diagnostics 
to compare the structures: the radial distribution function (rdf), the angular distribution 
function (adf) of nearest neighbors separated by a distance less than the first minimum of 
the rdf (usually r < 1.6 a), and bond order diagrams, which aid in identifying the global 
symmetry even if the symmetry elements are oriented in random directions. 

1. Radial Distribution Functions 

Figure ||a compares the radial distribution functions of the structures most commonly 
observed in the simulations. Typical of the amorphous structures is an asymmetric first 
peak, which appears to consist of two overlapping shells, followed by a second maximum in 
the range 1.7 < r < 2.7 a, which is the well-known double peak. The unusually sharp slope 
on the short- distance side indicates a well-defined second-nearest-neighbor distance, which 
is caused by the repulsive part of the maximum of the potential. The maximum at about 
1.9 a is formed by two opposite corners of a bipyramid consisting of two regular tetrahedra. 
A second remarkable sharp slope is found at the third maximum at about 2.7 o. Notice that 
the rdfs of the tcp structures obtained on cooling are completely indistinguishable from the 
amorphous case. 

The rdf of the cr-phase shows an example of the perfect square-triangle-rhombi-shield 
phases. The quasicrystalline and other crystalline and approximant phases differ only in the 
fine structure of the subpeaks. Evidently the tcp-phase rdf is a broadened envelope of the 
cr-phase rdf. 

The rdfs of the bcc and fee structures, however, are radically different. For bcc, the 
first maximum is indeed split, and in the second maximum the short- distance part is lower 



than the next peak. This maximum is now formed by the distances across the tetragonal 
octahedron in the bcc structure. At higher temperatures, when the peaks of the rdf are 
broadened and overlap, the rdf of the bcc-phase is similar to the rdf of the tcp-phases except 
that the weights of the two subpeaks of the second maximum between 1.7 a and 2.7 a are 
interchanged. 

The rdf of the fee crystal exhibits a peak at about 1.6 a. Since this is the position of 
the potential maximum, it is clear that the fec-phase is unstable at low pressures. Figure 
[5|b shows the transition from a bcc-phase to the fec-phase at k^T/e = 0.75. The small 
maximum at 1.5 a indicates formation of regular squares, which are characteristic of fee. 
The MD data are strikingly similar to theoretical predictions for the hard-sphere fee crystal 
(Fig.|c). 

The sequence of phases at low temperatures and increasing pressure becomes clear when 
we examine the rdfs: for bcc and the cx-phase, the first two overlapping atomic shells occupy 
the minimum of the potential. The next shell is beyond the maximum. For fee, the first 
shell is also in the minimum and the second shell is at the maximum. If the structures are 
compressed, the energies of bcc and the a-phase increase, since the second maximum of the 
rdf moves up the maximum of the potential. The energy of fee decreases, since the second 
maximum of the rdf moves down the potential maximum. 

2. Angular Distribution Functions 

The results for the angular distribution functions (Fig. |6|) are consistent with the results 
for the rdfs. The amorphous and tcp structures are indistinguishable. The tcp-phase adf 
is a broadened version of the adf of the a-phase. All of these phases, as well as the liquid 
show two maxima: a rather narrow extremum at small angles around 60° and a broad peak 
at about 120°. Both maxima indicate the existence of equilateral triangles. The adfs of bcc 
and fee are again completely different. Especially remarkable are the maxima at an angle of 
90°, indicating the existence of rectangles in these phases. For bcc these angles are formed 



by distances between atoms along the four-fold axis. 



3. Bond Order Diagrams 

We have seen that the amorphous and tcp structures obtained by cooling cannot be 
distinguished by the rdf and adf. Are they really different? As shown in Fig. [7], bond 
order diagrams can give a clear answer. The diagram for bcc (Fig. |7|a) is relatively simple, 
exhibiting seven maxima - four from the separation along the space diagonals and three from 
the distances along the four-fold directions. As the sample shown here was not perfect, thin 
bridges join the maxima. The liquid (not shown) is characterized by an isotropic distribution 
of the bonds covering the whole sphere homogeneously. The diagram for the tcp-phase 
(Fig. 0b) is somewhat more complicated, featuring an equator with twelve maxima indicating 
the presence of a quasicrystal. Other prominent features are two further circles of maxima 
at higher latitudes and two peaks at the poles. For non-perfect samples the distribution of 
the maxima is distorted, and the symmetry may not be dodecagonal. Occasional ring-like 
arrangements of overlapping maxima surrounded by further maxima suggest twinning and 
multi-grain samples. 

Comparing samples classified as amorphous or tcp-phases, we find that a continuous 
transition between the two may be possible. In Fig. |7|b the maxima can be seen quite 
clearly. In other samples, however, the maxima are almost obscured by a rather homogeneous 
background noise. Quenching and annealing improves the diagrams only marginally. It is 
possible that some of the amorphous samples actually consist of a number of micro-grains. 
If in fact the case, this would mean that the amorphous and tcp-phases have the same local 
arrangement of atoms, although the amorphous sample did not succeed in ordering globally. 

4- Real-Space Representation of the Structures 

Real space pictures (snapshots) of the samples also help to distinguish the structures. In 
most cases, the bcc samples look quite defect-free, with only the vacancies visible. Sometimes 



we find two differently oriented domains in the simulation box. Liquid samples obviously do 
not show any regularities. In the quenched amorphous structures, however, there are some- 
times partially ordered parts, underlining our claim that they contain micrograms. The tcp 
samples, as exemplified by Fig. 0c, have layered structures, which, if viewed perpendicu- 
lar to the layers, resemble the perfect sample in Fig. |3], characterized by centered ring-like 
structures formed by the 14- and 15-fold coordinated atoms. The quality of the pictures is 
often low, however, as the samples may be twinned or contain several grains. 

B. Ground-State Structures 

The equilibrium structure of a specific potential at a given temperature and pressure may 
be determined, in principle, by a global minimization of the Gibbs free energy G. In practice, 
however, this procedure is not feasible with MD simulations, since direct transitions between 
local minima are rarely observed. Even if special methods are used to switch between closely 
related structures like bcc and fee, there still may exist a free energy barrier high enough 
to prevent a transition [ |24}| . Instead of attempting to minimize G, one identifies promising 
structures, computes the thermodynamic functions, and compares them. An alternative is 
thermodynamic integration of P, starting at P = and integrating towards higher pressures. 

At T = 0, where entropy no longer affects stability, identifying the ground state simplifies 
considerably. Here the Helmholtz free energy F equals the internal energy U, which in turn 
equals the potential energy E pot , since the kinetic energy vanishes. Furthermore, the Gibbs 
free energy G becomes equal to the enthalpy H and the pressure P is determined by the 
virial equation, since the kinetic pressure k^T also vanishes. Common tangent constructions 
on curves of U(V,T — 0) vs. V yield the stability ranges of competing phases and curves of 
enthalpy H(P, T = 0) vs. P intersect at the phase coexistence pressures. 

We have calculated the ground-state energies by two independent methods. First, taking 
the perfect structures, we have computed lattice summations of the Dzugutov pair potential. 
Second, starting from a perfect structure, we have relaxed the system with the MD simulation 



program in an isothermal-isobaric ensemble, setting k-gT/e = 0.001 and P to the desired 
value. In this mode, molecular dynamics acts as a steepest-gradient optimization algorithm. 
The pressure is derived from the virial equation. In contrast to a lattice sum calculation, 
where only the volume is scaled, in the simulations all atoms may move independently. 
Therefore the results may (and do) differ slightly from the lattice sum calculations for 
perfect structures. 

Since the Dzugutov potential [Eq. ([!])] is isotropic and has a single minimum, it should 
favor densely packed structures if the volume is not restricted. An optimal packing in three 
dimensions would consist of regular tetrahedra, but such a packing does not exist. Now 
there are two choices to solve this dilemma: either introduce other coordination polyhedra, 
as in fee crystals, or use irregular tetrahedra, as in tcp phases. In the Frank-Kasper phases 
the coordination polyhedra are additionally restricted to deltahedra with five or six triangles 
meeting at a vertex. This condition is fulfilled for the icosahedron, and certain polyhedra 
with 14, 15 and 16 vertices. Although bec is tetrahedrally close-packed, it is not a Frank- 
Kasper phase since its coordination polyhedron (a rhombic dodecahedron) has vertices where 
only four triangles meet. 

In a first step towards identifying stable structures we study stacking variants, distort 
the phases mentioned above, and examine various tcp structures. The fee structure can 
be modified by stacking the densely packed layers differently. We find that the hexagonal 
close-packed (hep) and other stacking variants are considerably less stable than fee at high 
pressures where fee is more stable than bec and the a-phase. This result is remarkable since 
for the Lennard- Jones potential hep is known to be slightly more stable than fee ||25|| . 

Distortions of the bec phase along the principal symmetry axis always reduce the stability. 
The same happens for the a-phase if the layer distance is changed from the optimum at 
c/a = 1.03 (a is the edge length of the tiles and c the period along the z-axis). 

The Frank-Kasper phases are of two types: structures with 16-fold coordinated sites 
and structures without. The dodecagonal quasicrystal, its approximants, and crystalline 
variants are of the latter type, called the square-triangle class. Structures containing 16-fold 



(or higher) coordinated atoms have a lower stability. The sites with the high coordination 
numbers are too numerous and the potential energy increases because of strained bonds. 

We observe the same trend in the square-triangle class. The stability is lowest for the 
purely triangular Z-phase since the number of 15-fold sites is also considerable. The stability 
increases if the triangles are separated by squares, but is again rather low if the structure 
contains only squares as in the A-phase without 15-fold sites, perhaps because the A-phase 
has full cubic symmetry and is therefore more rigid than the other structures. The a-phase, 
on the other hand, is more stable than the H-phase, since it contains only pairs of triangles 
instead of rows. More complicated crystalline phases, approximants, and the quasicrystals 
all contain mixtures of squares and triangles in different arrangements. These structures are 
all inferior to the a-phase since they must contain larger conglomerates of triangles. 

In a second step towards identifying ground-state structures, we survey the published 



crystallographic structures. From the lists in Refs. [26-28 1 , a variety of structures have been 
selected according to the following criteria: 

1. Coordination numbers between 10 and 15. 

2. Derivatives of tcp structures. 

3. Derivatives of bcc, especially vacancy-ordered structures. 

4. Quasicrystal approximants. 

5. Icosahedral coordination shells. 

For each structure examined the required crystallographic data were taken from Ref. || . A 
full list of the structures is given in the Appendix. 

Assembling the results, the following picture emerges at T = (Fig. |8|): the bcc-phase has 
the absolute minimum potential energy at a density of pa 3 = 0.866. The a-phase acts as the 
lower bound for all the square-triangle phases, being minimal at pa 3 = 0.879. The fcc-phase 
has a potential energy minimum at pa 3 = 1.013. A common tangent construction shows that 



bcc is stable up to pa 3 = 0.887, and fee is stable above per 3 = 1.057. The relaxed cr-phase is 
stable, as determined by MD, only within the narrow interval 0.887 < per 3 < 1.057, whereas 
the ideal cr-phase is never stable, according to a simple lattice sum calculation. Intersections 
of the enthalpy curves yield the stability ranges in terms of pressure. The most stable 
structures are bcc for Pa 3 /e < 1.70, cr for 1.70 < Pa 3 /e < 2.85, and fee for Pa 3 /e > 2.85. 
The sequence of ground-state structures with increasing pressure (and density) is therefore: 
bcc - a - fee. The properties of the various structures are summarized in Table |. 

Vacancy-ordered phases are more stable than pure bcc at densities down to pa 3 = 0.7 
(Fig. ||b). In the range 0.6 < pa 3 < 0.7 the lowest potential energy is attained by a disordered 
phase formed upon annealing the NiTi2 approximant phase. 



C. Finite- Temperature Phases 

From the ground-state calculations we have determined the stable structures at T = 0. 
Insight into the topology of the phase diagram at finite temperatures now can be obtained 
by observing phase transitions between the melt and the solid upon heating, cooling, and 
compression. As noted in Sec. [VTJ, it is not possible to determine the relative thermodynamic 
stabilities of two competing phases directly from conventional MD simulations because of the 
difficulty of computing the entropic contribution to the free energy. However, by observing 
the temperature and pressure at which a phase becomes unstable, it is possible to establish 
limits of mechanical stability. The results of our MD stability analysis are consolidated in 
Fig. Ilk. 



1. Heating Simulations 

If the bcc and cr-phase solids are heated at low pressure, the energy and enthalpy for 
bcc remain always lower than those for the cr-phase. At higher pressures, the energy of 
the cr-phase drops below that of bcc. The differences between the enthalpies at higher 



temperatures, however, are smaller than their fluctuations, such that the relative stabilities 
of bcc and the cx-phase cannot be resolved. 



The determination of the melting line has been discussed in detail in Ref. [29]. Here we 
present only a brief summary. The phase transition line was determined by preparing a solid 
at ksT/e = or 0.4 at fixed pressure and heating it continuously at rates of k^ST /e=0.001 or 
0.002 per timestep until melting was observed. The criterion for melting was the divergence 
of the mean-square displacement. At the same temperature a sudden rise in the potential 
energy and an associated drop in the density were observed. Similar simulations have been 
carried out at constant volume starting at k B T/e = 0.001 and Pa 3 /e = 0.001. We emphasize 
that the transition line thus obtained is not strictly the equilibrium melting line, since 
with periodic boundary conditions the sample has no surface at which melting could start 
and a two-phase coexistence is not possible because the samples are too small. The solid- 
fluid transition lines differ only slightly for bcc and the cx-phase. The fee crystal melts at 
somewhat higher temperatures. At high pressures, however, the cx-phase becomes unstable 
at considerably lower temperatures than bcc and fee. 

The transition to the fluid also can be determined by expanding a solid at constant tem- 
perature starting from high densities. The transition line obtained in this way for bcc crystals 
is the same as that determined by heating within the statistical fluctuations (Fig. |i~0|a). As 
noted in Ref. ||29|| , only one fluid phase is observed and no transition between a liquid and 
a vapor phase could be found. 



2. Cooling Simulations 

Cooling simulations were carried out in a manner similar to the heating runs. Starting 
samples were obtained from solids equilibrated at high temperatures. The cooling rate was 
kv5T/e = 0.002 per time step. Similar to melting, the freezing transition is delayed, now 
because critical nuclei first must be formed, and subsequent large-scale reordering of atoms 
may be necessary. 



We find that the phase nucleating at pressures above Per 3 / e = 5 always has bec symmetry. 
If the temperature is lowered to about k-gT/e = 0.7 — 1.0, bec becomes unstable relative to 
fee at pressures above Pa 3 /e = 20. Although a complete transition to fee cannot be achieved 
with our simulation method, we observe a clear indication that fee is the preferred structure. 
In the radial distribution function (Fig. |5]a) we observe the emergence and growth of a new 
peak between the first and the second peaks at about v2 times the nearest-neighbor distance, 
which signals the formation of regular squares in the close-packed crystal structures. 

Below Pcr 3 /e = 5 we do not observe a typical freezing transition with a jump in potential 
energy, but only a sharp kink, reminiscent of a glass transition. The nucleating structures 
are partially ordered and possess features typical of the tcp structure, namely layering, ring- 
like structures, and the Frank-Kasper polyhedra (see Fig. |7]c). In the following, we refer to 
such structures as the tcp-phase. Although sometimes dodecagonal, the tcp structures often 
do not have a perfect symmetry, and thus may have varying degrees of crystallinity. 

Because of the maximum in the potential, it was not possible in general to obtain perfect 
samples. If the pressure is too low, there is insufficient cohesion to compactify the samples. 
This is clearly seen by comparing runs with different cooling rates. However, if equilibrated 
for a longer time, the samples eventually become much denser. 

The density ranges for stability have been obtained by cooling at constant volume. For 
the 500-atom sample we obtain the boundary between the formation of the bec and the 
tcp-phase at pa 3 = 0.87, independent of the cooling rate up to k^ST/e = 0.0005 per time 
step. For the 1024-atom sample, however, the boundary is shifted to per 3 = 0.84 and is 
observed at the first time for k^ST/e = 0.00025 per time step. This is remarkable, since the 
minimum of the a-phase lies at about pa 3 = 0.9. 

The formation of the crystalline structures also depends on the sample sizes and the 
cooling rates. A sample with 250 atoms and a constant density of pa 3 = 0.865 froze to bec 
at a cooling rate of ksST/e = 0.001. For 500 atoms we had to reduce the cooling rate by a 
factor of one half, and for 1024 atoms a cooling rate of k-gST/e = 0.00025 per time step was 
necessary to obtain a perfect bec phase, although partial bec ordering was already observed 



at twice this rate. 

It is easier to obtain the tcp-phase in a constant volume simulation (as Dzugutov did) 
rather than in a constant pressure ensemble. To some extent, the nucleated structures can be 
annealed also at constant volume. However, most of the defects, especially different domains, 
cannot be so removed. Annealing at constant pressure also turns out to be ineffectual. 

The transition from fluid to solid also may be observed by compressing the fluid at a 
constant pressure gradient of 5Pa 3 /e = ±0.1. The transition curve is the same as for cooling 
(Fig. |lO|a), the collapsed structures being again bcc, at least for k B T/e = 0.6, 0.8, 1.0, 1.5, 
and 2.5. 

Between the melting and freezing curves we observe a broad hysteresis region, within 
which the thermodynamic phase transition should occur. The reason for the broad hysteresis 
region is the peculiar form of the potential. The maximum strongly inhibits freezing and 



collapsing of the structure, indeed as intended by Dzugutov [[30 



The structures generated by cooling the samples contain free volumes even if the density 
or pressure during nucleation is high. In the case of constant volume cooling the reason 
is obvious, since the volume of the nucleating regions shrinks with temperature. However, 
constant pressure cooling also generates free volumes, even at high pressures, since at the 
onset of nucleation the frozen domains have a higher density than the liquid. The rigidity 
of the solid prevents the simulation box from contracting fast enough. With the methods 



described in Sec. [II Q , we can show that the free volumes in the ordered phases are mostly 
vacancies. 

If the simulation samples are quenched to T = and P = 0, and the vacancies are 
filled with atoms, we find that the density of bcc rises to per 3 = 0.864 ± 0.005, whereas the 
densities of the tcp-phase remain at about pa 3 = 0.847 ± 0.005. Although the densities of 
the bcc samples are close to the ideal value of 0.8638, the densities of the tcp-phase are far 
lower than the ideal value at the potential energy minimum (per 3 = 0.881). 



3. Compression Simulations 



If the structures are compressed at fixed temperature, bcc destabilizes first. One might 
therefore expect bcc to be stable only at relatively low temperatures. However, this would 
contradict the cooling simulations (Sec. |V C 2| ) , which yield a bcc structure. A full picture can 
be obtained only by calculating the Gibbs free energy, since it may be kinetically favorable 
for the system to nucleate bcc crystallites. 

At high pressures the stable structure is clearly fee, which has the lowest energy and 
enthalpy. Upon compression, this close-packed structure remains stable, and radial distri- 
bution functions of the decaying bcc and a-phase structures show new peaks characteristic 
of fee. 



D. Theoretical Predictions 

For comparison with the MD simulation data, we have applied the perturbation theory 
described in Sec. |T| to predict the thermodynamic phase behavior of the Dzugutov-potential 
system. For the fluid phase and selected solid structures, free energies were calculated and 
a coexistence analysis performed. Our choice of structures was dictated by the structures 
actually observed in the simulations. Figure ^| compares the HS part of the free energy for the 
fee, bcc, and cx-phase structures. Also shown, for comparison, are corresponding Monte Carlo 



simulation data from Ref. ||19|| . From the maximum HS volume fractions of these structures - 
respectively 0.74, 0.68, and 0.53 - stability of the HS solid is seen to be strongly influenced by 
packing efficiency. Therefore, at high temperatures and pressures, where entropy dominates 
the free energy and the system behaves as an effective HS system, the structures that 
are more efficiently packed are favored. As temperature and pressure decrease, internal 
energy makes an increasing contributon to the free energy. As illustrated by the neighbor 
distance histograms (Fig. P and the corresponding hard-sphere rdfs (Fig. |5|c), the first 
few coordination shells of the bcc and cr-phase structures are more commensurate with the 



attractive part of the Dzugutov potential than those of the fee crystal, favoring these more 
loosely packed structures over close-packed fee. 

Constructing Maxwell common tangents to curves of free energy per volume vs. density, 
thus ensuring equality of chemical potentials and pressures in coexisting phases, we have 
mapped out the phase diagram of the system. Projections onto the P — T and T — p 



planes are shown in Figs. |I0| b and [II], respectively. As anticipated, the stable solid at high 
pressures is the fee crystal, while the bee crystal is only metastable relative to fee (long- 
dashed curve in Fig. |lO|b). Aside from fee, bee, and a-phase, we have also considered several 
tcp structures observed in the simulations and rational approximants to layered dodecagonal 
quasicrystals. The tcp and quasicrystal structures, however, were found to be at best only 
metastable relative to the crystal structures. At T = 0, lattice-sum calculations of ground- 
state energies (Fig. ||a) show that bee is the stable structure for Pa 3 /e < 2.66. From this 
known limit, we postulate that bee is also the stable solid structure at low P for small but 
finite temperatures. The perturbation theory being of uncertain accuracy for k^T je < 0.5, 
we further postulate an extrapolation of the fluid-fee phase boundary to zero pressure. This 
confines the stable bee phase to a small pocket in the lower-left corner of the P — T diagram. 



VI. DISCUSSION AND CONCLUSIONS 

The pressure-temperature phase diagram of the Dzugutov potential obtained by MD 
simulation is surprisingly rich. At low pressures and temperatures, the bec-phase is stable, 
followed, with increasing pressure, by the cx-phase and by fee (Fig. p!0|a). The bee crystal 
is nucleated from the fluid for sufficiently slow cooling rates and sufficiently high density or 
pressure. It is also obtained by compressing the fluid. Below Pa 3 /e = 5 or pa 3 = 0.85, tcp 
structures, including the dodecagonal quasicrystal, are formed. The cooling scenario may 
be summarized as follows: 

1. At high cooling rates a glass is formed, which may be transformed into tcp or bec 
solids by annealing at ksT/e = 0.4 — 0.5. 



2. At lower rates the fluid has sufficient time to reorder locally and crystallizes into a 
bcc crystal. The bcc structure being relatively simple, the samples are in most cases 
perfect except for vacancies. 

3. At sufficiently low density or pressure, a tcp structure is generated. Characterized 
by layered symmetry, the tcp-phase can be formed for a wide variety of coordination 
polyhedra and many energetically similar configurations. 

4. If the bcc crystal is cooled at high pressures or compressed at low temperatures, it 
transforms, at least partially, into the fee structure. 

The stability of the lowest energy tcp-phase, namely the cr-phase, relative to bcc, could 
not be determined precisely by our simulations. Heating at low pressure shows that the 
energy and enthalpy of bcc are always lower than those of the cr-phase. Comparing the 
two phases at higher pressures and temperatures shows that the difference in enthalpy is no 
longer significant, but the energy of the cr-phase becomes lower than that of bcc at higher 
pressures and temperatures. Furthermore, upon compression bcc becomes unstable at lower 
pressure than the cr-phase. The detailed topology of the phase diagram in Fig. |TD|a is still 
not completely clear and further simulations are necessary to compute the phase boundaries 
exactly. At low temperatures bcc appears at lower pressures than the cr-phase, whereas in 
the cooling simulations bcc is formed at higher pressures compared to the tcp structures. It 
may be, as speculated by Dzugutov [|J, that entropy lowers the free energy of the cr-phase, 
and especially of the quasicrystal, thus leading to a stable tcp or quasicrystalline state at 
higher temperatures. 

Our theoretical calculations for a selection of perfect solid structures suggest that the 
thermodynamically stable solid phases of the Dzugutov-potential system are limited to fee 
and bcc crystals. Lattice-sum calculations at T = show that the cr-phase is almost de- 
generate with, though of slightly higher energy than, the bcc crystal. At high temperatures 
(ksT ^> e), where the attractive well in the potential plays only a minor role, packing effi- 
ciency strongly disfavors the cr-phase relative to both fee and bcc crystals. At intermediate 



temperatures (k^T ~ e), perturbation theory predicts fee and bee crystals to be always 
more stable than the relatively loosely packed a-phase. Thus the a-phase appears nowhere 
in the P — T phase diagram and the bee crystal appears only at low P and T. 

It may be, of course, that the first-order perturbation theory lacks sufficient accuracy 
to conclusively resolve relative stabilities of such closely competing phases. In particular, 
fluctuations in the total perturbation energy, being stronger for a disorder fluid than for 
ordered solids, render the theory inherently less acurate for the fluid phase. Moreover, the 
mean-field neglect of next-nearest-neighbor correlations in the HS rdf is less justifiable for 
more open structures, such as bee and the a-phase, than it is for the close-packed fee crystal. 
Such correlations, if significant, would tend to lower the free energies of the open structures 
and might possibly influence the order of stabilities. 

Nevertheless, it should be emphasized that the predicted phase diagram is not neces- 
sarily at odds with the simulation data, if the tcp and high-T bee phases, observed in the 
simulations, are regarded as metastable with respect to fee. Conceivably, for kinetic rea- 
sons, the supercooled fluid first nucleates a metastable bee crystallite, which upon growth 
transforms into a stable fee crystal. In fact, such behavior has been predicted from a 
general density- wave instability argument |nj, and has been observed in simulations of a 
supercooled Lennard- Jones fluid [fffil - Furthermore, whereas the theory has been applied to 
perfect structures, the simulations often result in solids replete with defects. We may conjec- 
ture, therefore, that the defects in the a- and tcp-phases observed in the simulations serve 
to improve the packing efficiencies {e.g., by increasing nearest-neighbor distances), while 
approximately preserving the average coordination distances, thereby confering energetic 
advantage over the fee structure. 

After examining a wide variety of solid structures as candidates for stable phases of 
the Dzugutov potential, we are drawn to conclude that only such simple structures as bec 
or fee are competitive. A possible exception is the a-phase, a tetrahedrally close-packed 
structure, although one of the simplest examples of its class. These results appear to place 
the Dzugutov potential in line with the Yukawa, generalized Lennard- Jones, Rubidium, 



and Morse potentials, all of which favor bcc, fee, or hep crystals, supporting the general 
principle that simple pair potentials tend to favor simple structures. Nevertheless, it remains 
conceivable that quasicrystals and other complex structures, such as A15, Z, or H, might 
gain stability through modifications of the Dzugutov potential. Indeed previous work has 
identified somewhat related pair potentials for simple metals - albeit with no counterpart in 
the periodic table - for which stable icosahedral quasicrystals have been predicted [33,34 



Future work along these lines could examine variations of the Dzugutov potential in attempts 
to modify the relative stabilities of the tcp-phases and dodecagonal quasicrystals. 
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APPENDIX 

The following lists contain the structures investigated as possible ground states in the 
MD simulations. The notation is as found in Ref. ||. 

Stacking variants of fcc(abc): hcp(ab), abcacb, abebeb. 

Frank-Kasper phases: A15 or /3-W (cP8 Ch^Si), Z (hp7 Zr 3 Al 4 ), cr-phase or /3-U (tP30 
Cr 46 Fe 54 ), H, J, F, K', C15 (cF24 MgCu 2 ), T (cI162 Mg 32 (ZnAl) 4 9), n (hR13 Fe 7 W 6 ), D 
(V 2 6Fe44Si 3 o). 

Square-triangle phases: A, H, Z, cr-phase, J, F, K', inflations of A and Z, doubling of 
the tiles of cr-phase, a phase with a and h sites mixed, tetragonal approximants of the 
quasicrystal with 23, 36, 172 and 836 tiles. 

Vacancy ordered phases derived from bcc: cPl, cIlO /3-Hg 4 Pt, cF12 CaF 2 , tC14 AsPdsTl, 
cI52 Cu 5 Zn 8 , cF120 Pd 4 _ x Te, cF120 Scnlr 4 , cF116 Mn 23 Th 6 , cF88 Bi 4 Cu 4 Mn 3 , cI112 



Ga 4 Ni 3 . 

The other structures are: oP12 Co 2 Si, 0PI6 AlDy, 0C8 BCr, oCIO AlFe 2 B 2 , oC12 Ge 2 Th, 
oC12 Si 2 Zr, 0CI6 BCMo 2 , 0CI6 Ga 3 Pt 5 , 0CI6 HgNa, oC28 Al 6 Mn, oIlO B 2 CoW 2 , oI12 
Gd 2 Si 3 , 0II6 BMo, oI20 AI4U, oF24 Si 2 Ti, oF48 CuMg 2 , tpl4 Hg 5 Mn 2 , tP20 Al 2 Gd 3 , tP30 
AlNb 2 , tI12 Al 2 Cu, tI12 Si 2 Th, tI16 BMo, tI28 MnU 6 , tI32 Si 3 W 5 , hP3 Cd 2 Ce, hP3 A1B 2 , 
hP5 Al 3 Ni 2 , hP6 InNi 2 , hP6 Caln 2 , hR7 B 5 Mo 2 , cP8 FeSi, cP20 Mn, cP39 Mg 2 Znn, cP138 
Al 9 Mn 2 Sin, cP140 (AlSi) 58 Mni 2 , cI12 Ga, cI26 A1 12 W, cI58 Mn, cI76 Cui 5 Si 4 , cF96 NiTi 2 . 
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TABLES 



TABLE I. Ground-state structure and reduced density pa 3 at which the potential energy per 
atom E pot is minimal. The upper half contains MD results at k^T/e = 0.001, the lower half lattice 
sum calculations for perfect structures. 
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1.013 


-2.19 
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1.057 


-2.22 


bee 
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FIGURES 




FIG. 1. Dzugutov pair potential together with histograms of neighbor distances for fee (black 
bars), bec (gray bars), and u-phase (unshaded bars) crystals at reduced density per 3 = 0.9. Heights 
of bars are proportional to numbers of neighbors. 




FIG. 2. The basic tiles of the dodecagonal model: square, triangle, rhombus, shield, twofold 
symmetric hexagon. The dotted atoms are placed in j4-layers z = 1/4 and 3/4, the white atoms in 
-Bdayers at z = and the black atoms at Z = 1/2. All tiles can also occur with black and white 
atoms exchanged, depending on their orientation. The twofold symmetric hexagon is unstable and 
does not occur in our tilings. 



FIG. 3. Patch of a decorated square-triangle quasicrystal. The dotted atoms are at z = 1/4 
and 3/4 in the A layer, the white atoms are at z = in the B layer and the black atoms at z = 1/2 
in the B layer. 




FIG. 4. Vertex configurations of the crystalline phases with a single vertex configuration. From 
left to right: A-, Z-, H- and <7-phase or a, z, h, and a vertex, respectively. 
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FIG. 5. Radial distribution functions for various systems, (a) From bottom to top: fee, bee, 
cr-phase, tcp-phase, amorphous phase. The samples have been expanded to Pa 3 /e = 0.001 and 
quenched to T = 0. Vertical scale is in arbitrary units, all curves being scaled to the same 
maximum, (b) A sample obtained from cooling at Pcr 3 /e = 25 and k^T/e = 0.75. The shoulder 
at r = 1.5 a indicates the onset of a transition to fee; (c) Hard-sphere solid at reduced density 
pa 3 = 1.0, computed from Eqs. (|i~3|)-(|l5|): fee crystal (solid curve), bec crystal (long-dashed curve), 
and cr-phase (short-dashed curve). 
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FIG. 6. Angular distribution functions. From bottom to top: fee, bee, cr-phase, tcp-phase, 
amorphous phase. The samples have been expanded to Pa 3 /e = 0.001 and quenched to T = 0. 
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FIG. 7. (a) Spherical projection of the nearest-neighbor vectors of a non-perfect bcc sample 
generated by cooling. Three views of the sphere along perpendicular directions are given; (b) 
Spherical projections of the nearest-neighbor vectors of a tcp sample obtained by cooling. Three 
perpendicular views are given. The poles of the sphere are marked by arrows, and the equator is 
indicated with dashes. The six maxima along the equator represent the twelve-fold symmetry; (c) 
Projection of a tcp sample obtained by cooling at Pa 3 /e = 3.5 and k%,T/e = 0.55. The sample has 
been expanded to P<r 3 /e = 0.001 and quenched to T = after nucleation. 




FIG. 8. Ground-state energy per unit volume vs. density, in reduced units, (a) Results of 
lattice summation of the Dzugutov pair potential for ideal fee crystal (solid curve), bee crystal 
(long-dashed curve), and cr-phase (short-dashed curve); (b) Results of MD simulation, with struc- 
tural relaxation, for cr-phase (solid curve), fee (long-dashed curve), and bec (short-dashed curve); 
(c) MD simulation data for the bec vacancy phases and the low density amorphous structures. The 
lowest minimum at the right is bec followed by CusZns, Pd4_ x Te, Mn23Tli6, and Ga4Ni3. The 
next minimum belongs to the amorphous phase formed from NiTi2. The remaining double-dashed 
curve is the amorphous phase generated by cooling the melt. 
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FIG. 9. Free energy per volume vs. density, in reduced units, for the reference hard-sphere 
solid, computed from Eqs. (^), (||), and (|T2|). Curves have the same meaning as in Fig. ||a. Circular 



and square symbols are Monte Carlo simulation data, from Ref. [19], for hard-sphere bec and fee 
crystals, respectively. 
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FIG. 10. Pressure-temperature phase diagram, (a) MD simulation results: The instability lines 
denote boundaries where respective structures are destabilized if compressed to high pressures. 
Capital letters mark phases formed by cooling simulations, lowercase letters phases obtained by 
ground-state structure calculations. Crosses characterize region where tcp-phase is found in cooling 
simulations. The region between the melting/expansion transition line and the cooling/compression 
transition line is the hysteresis region; (b) Perturbation theory predictions: Phase boundaries are 
shown between fluid and fee crystal (solid curve) and between fluid and metastable bec crystal 
(long-dashed curve). Short-dashed curves are postulated extrapolations to low P and T. 
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FIG. 11. Predictions of perturbation theory for the fluid-solid phase diagram in the tempera- 
ture-density (T — p) plane. For k^T/e > 0.5, theory predicts the fee crystal to be the only stable 
solid phase. 



